#include <math.h>
#ifndef FFTCONST_H
#include "fftconst.h"
#endif // FFTCONST_H
//此代码来源《数字信号处理C语言程序集》殷福亮、宋爱军，沈阳：辽宁科学技术出版社，1997.7
//数组x存储时域序列的实部，数组y存储时域序列的虚部
//n代表n点FFT，sign=0为FFT，sign=1为IFFT

void FFT (myint x[], myint y[], int n, int sign)
{
    int i, j, k, l, m, n1, n2;
    double c, c1, e, s, s1, t, tr, ti;

//unsigned char c,c1,e,s,s1,t,tr,ti;
//unsigned char i,j,k,l,m,n1,n2;
    /*----计算分解的级数m=log2(n)----*/
    for (j = 1, i = 1; i < n; i++)
    {
        m = i;
        j = 2 * j;
        if (j >= n)
            break;
    }
    /*----按照倒位序重新排列原信号----*/
//即：计算蝶形图的输入下标（码位倒读）
    n1 = n - 1;
    for (j = 0, i = 0; i < n1; i++)
    {
        if (i < j)
        {
            tr = x[j];
            ti = y[j];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = tr;
            y[i] = ti;
        }
        k = n / 2;
        while (k < (j + 1))
        {
            j = j - k;
            k = k / 2;
        }
        j = j + k;
    }

//计算每一级的输出，l为某一级，i为同一级的不同群，使用同一内存（即位运算）
    n1 = 1;
    for (l =0; l < m; l++)
    {
        n1 = 2 * n1;//n1=2*n1=2^m代表第m级每个分组所含节点数,m=log2(n)
        n2 = n1 / 2;//n2代表第m级每个分组所含碟形单元数,//同时它也表示每个碟形单元上下节点之间的距离
        e = PI / n2;
        c = 1.0;
        s = 0.0;
        c1 = cos (e);
        s1 = sin (e);
        if (0 == sign)
            s1 = -s1;
        /*----碟形运算----*/
        for (j = 0; j < n2; j++)
        {
            //遍历每个分组，分组总数为n/n1
            for (i = j; i < n; i += n1)
            {
                k = i + n2;//i,k分别代表一个碟形单元的上、下节点编号
                tr = c * x[k] - s * y[k];//wn=Wnr
                ti = c * y[k] + s * x[k];
                x[k] = x[i] - tr;
                y[k] = y[i] - ti;
                x[i] = x[i] + tr;
                y[i] = y[i] + ti;
            }
            t = c;
            c = c * c1 - s * s1;
            s = t * s1 + s * c1;
        }
    }
//如果求傅里叶逆变换IFFT，再除以n
    if (sign == 1)
    {
        for (i = 0; i < n; i++)
        {
            x[i] /= n;
            y[i] /= n;
        }
    }
}

void
IFFT (myint x[], myint y[], int n, int sign)
{
    fix_FFT(&x,&y,n,sign);
}


void fix_FFT (myint x[], myint y[], int n, int sign)
{//ligb20170423 这是我自己根据FFT函数修改的定点整数FFT，因为从网络获得的开源代码不能计算幅角
    int i, j, k, l, m, n1, n2,n3;
    double c, c1, e, s, s1, t;
    double tr, ti;
    char k1;//记录
    double c2,s2;
    n3=1<<((sizeof(Sinewave[0])<<3)-1);

//unsigned char c,c1,e,s,s1,t,tr,ti;
//unsigned char i,j,k,l,m,n1,n2;
    /*----计算分解的级数m=log2(n)----*/
    for (j = 1, i = 1; i < n; i++)
    {
        m = i;
        j <<=1;
        if (j >= n)
            break;
    }
    /*----按照倒位序重新排列原信号----*/
//即：计算蝶形图的输入下标（码位倒读）
    n1 = n - 1;
    for (j = 0, i = 0; i < n1; i++)
    {
        if (i < j)
        {
            tr = x[j];
            ti = y[j];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = tr;
            y[i] = ti;
        }
        k = n >>1;
        while (k < (j + 1))
        {
            j = j - k;
            k >>= 1;
        }
        j = j + k;
    }

//计算每一级的输出，l为某一级，i为同一级的不同群，使用同一内存（即位运算）
    n1 = 1;
    for (l = 0; l < m; l++)
    {
         n2 = n1;//n2代表第m级每个分组所含碟形单元数,//同时它也表示每个碟形单元上下节点之间的距离
        n1 <<= 1;//n1=2*n1=2^m代表第m级每个分组所含节点数,m=log2(n)
        e = PI / n2;
        c = 1.0;
        s = 0.0;
        k1=N_WAVE>>(l+1);
        c1=Sinewave[k1+N_WAVE/4-1]/(double)n3;
        s1=Sinewave[k1-1]/(double)n3;
        c2=c1;
        s2=s1;
        c1 = cos (e);
        s1 = sin (e);
        if (0 == sign)
            s1 = -s1;
        /*----碟形运算----*/
        for (j = 0; j < n2; j++)
        {
            //遍历每个分组，分组总数为n/n1
            for (i = j; i < n; i += n1)
            {
                k = i + n2;//i,k分别代表一个碟形单元的上、下节点编号
                tr = c * x[k] - s * y[k];//wn=Wnr
                ti = c * y[k] + s * x[k];
                tr;//因为旋转因子存储的是扩大了32767倍的整数值，所以要缩小
                ti;
                x[k] = x[i] - tr;
                y[k] = y[i] - ti;
                x[i] = x[i] + tr;
                y[i] = y[i] + ti;
            }
            t = c;
            c = c * c1 - s * s1;
            s = t * s1 + s * c1;
        }
    }
//如果求傅里叶逆变换IFFT，再除以n
    if (sign == 1)
    {
        for (i = 0; i < n; i++)
        {
            x[i] /= n;
            y[i] /= n;
        }
    }
}
